% run this file to recreate Figure 7 ("Impacts of Permanent
% Shocks : full Information versus learning")

% The figure includes the IRFs for 6 different cases: i) full information,
% ii) learning - early 2000s and iii) learning - late 2000s both for for Oil-supply
% shocks and TFP (neutral technology) shocks

comp = [0 1];
track = 0;

% IRFs for full information
% Calibration of oil supply shocks, early 2000s
rhoe_p = 0.99999;   % serial correlation of persistent energy shocks
sige_p =  0.0655*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.8158;   % serial correlation of transitory energy shocks
sige_t = 0.01;    % std. deviation of innovations to those shocks

%Calibration for TFP shocks "early 2000s"
rhoz_p = 0.99999;   % serial correlation of persistent productivity shocks
sigz_p = 0.2707*0.3694*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.8466;   % serial correlation of transitory productivity shocks
sigz_t = 0.3694*0.01;    % std. deviation of innovations to those shocks

whatinfo = 1;
% whatinfo = 1 implies "complete information" solution as defined in the paper; whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present  
answlearn = 1;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

gainkeep = []; %this is not used here but in other programs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1- Compute the "complete information" IRFs.

% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);


% compute 40-periods IRFs for all variables following a permanent oil-supply shock
nM = rows(M);
nPI = rows(PI);
S = zeros(1,nM)';
S(6) = -0.01271; % Shock magnitude such that the full-information impact response of oil prices is 10%
nir = 50;
IRB = zeros(nir,nPI+1);




for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR1 = IRB; % save responses

%now compute responses for shock to TFP
S = zeros(1,nM)';
S(8) = 0.0390; % Shock magnitude such that the full-information impact response of oil prices is 10%
nir = 50;
IRB = zeros(nir,nPI+1);

for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR2 = IRB; % save responses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2- Compute the "learning - early 2000s" IRFs.

whatinfo = 2;
param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];


lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

% compute IRFs for all variables following a permanent oil-supply shock
nM = rows(M);
nPI = rows(PI);
S = zeros(1,nM)';
S(6) = -0.01271;
nir = 50;
IRB = zeros(nir,nPI+1);

for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR3 = IRB; % save responses

%now compute responses for the shocks to TFP
S = zeros(1,nM)';
S(8) = 0.0390;
nir = 50;
IRB = zeros(nir,nPI+1);

for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR4 = IRB; % save responses 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3- Compute the "learning - late 2000s" IRFs.


% change calibration of shocks to the "late 2000s" values
rhoe_p = 0.99999;   % serial correlation of persistent energy shocks
sige_p = 0.2408*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.8294;   % serial correlation of transitory energy shocks
sige_t = 0.01;    % std. deviation of innovations to those shocks
% "late 2000s"
rhoz_p = 0.99999;   % serial correlation of persistent TFP shocks
sigz_p = 2.4926*0.2592*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.8463;   % serial correlation of transitory TFP shocks
sigz_t = 0.2592*0.01;    % std. deviation of innovations to those shocks

whatinfo = 2;

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];


% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

% compute 40 periods IRFs for all variables following a permanent oil-supply shock
nM = rows(M);
nPI = rows(PI);
S = zeros(1,nM)';
S(6) = -0.01271;
nir = 50;
IRB = zeros(nir,nPI+1);

for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR5 = IRB; % % save responses 

%now compute responses following shocks to TFP
S = zeros(1,nM)';
S(8) = 0.0390;
nir = 50;
IRB = zeros(nir,nPI+1);

for i=1:nir
   %Index of Period;
   IRB(i,1) = i-1;
   %Observables
   IRB(i,2:nPI+1)=(PI*S)';
   % Update the state vector:;
   S = M*S;  
end
IR6 = IRB; % save responses 

%Express all IRF is percentage points for graphiing purposes
IR1 = 100*IR1;
IR2 = 100*IR2;
IR3 = 100*IR3;
IR4 = 100*IR4;
IR5 = 100*IR5;
IR6 = 100*IR6;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Now ready to print responses

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(7)
% a 4 by 2 graph with oil-supply (left panel) and TPF (right panel) shocks side
% by side

trait1 = 'b-';
trait2 = 'r:';
trait3 = 'gx';

legende1 = 'Full Information';
legende2 = 'Learning, Early 2000s';
legende3 = 'Learning, Late 2000s';

myylabelshort = 'Dev. from s.s.  ';
mysize = 15;
myxlabel = 'Quarters';

fmt = 5.3;
epais = 4.5;


TT1 = 40;
dt1 = (0:TT1-1)';

subplot(421),plot(dt1,IR1(1:TT1,1+iy),trait1,dt1,IR3(1:TT1,1+iy),trait2,dt1,IR5(1:TT1,1+iy),trait3,'LineWidth',epais),
title('GDP','FontWeight','bold','FontSize',mysize)
%axis([0 TT2 -20 -15])
%xlabel(myxlabel,'FontSize',mysize3,'FontWeight','bold')
ylabel(myylabelshort,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);

h1 = text(12.0,-0.05,'Oil Supply Shocks');
h1.Interpreter = 'latex';
h1.FontWeight = 'bold';
h1.FontAngle = 'italic';
h1.FontSize = 20;


subplot(423),plot(dt1,IR1(1:TT1,1+is),trait1,dt1,IR3(1:TT1,1+is),trait2,dt1,IR5(1:TT1,1+is),trait3,'LineWidth',epais),
%axis([0 TT2 -1 1])
title('Oil Inventories','FontWeight','bold','FontSize',mysize)
%xlabel(myxlabel,'FontSize',mysize3,'FontWeight','bold')
ylabel(myylabelshort,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);

subplot(425),plot(dt1,IR1(1:TT1,1+ipo),trait1,dt1,IR3(1:TT1,1+ipo),trait2,dt1,IR5(1:TT1,1+ipo),trait3,'LineWidth',epais),
title('Oil Prices','FontWeight','bold','FontSize',mysize)
%axis([0 TT2 -55 -30])
%xlabel(myxlabel,'FontSize',mysize,'FontWeight','bold')
ylabel(myylabelshort,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);

subplot(427),plot(dt1,IR1(1:TT1,1+ipo)-IR1(1:TT1,1+ipof4),trait1,dt1,IR3(1:TT1,1+ipo)-IR3(1:TT1,1+ipof4),trait2,dt1,IR5(1:TT1,1+ipo)-IR5(1:TT1,1+ipof4),trait3,'LineWidth',epais),
title('Spot minus Futures','FontWeight','bold','FontSize',mysize)
%axis([0 TT2 -0.5 1.0])
xlabel(myxlabel,'FontSize',mysize,'FontWeight','bold')
ylabel(myylabelshort,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);



subplot(422),plot(dt1,IR2(1:TT1,1+iy),trait1,dt1,IR4(1:TT1,1+iy),trait2,dt1,IR6(1:TT1,1+iy),trait3,'LineWidth',epais),
title('GDP','FontWeight','bold','FontSize',mysize)
%axis([0 TT2 -20 -15])
%xlabel(myxlabel,'FontSize',mysize3,'FontWeight','bold')
ylabel(myylabelshort,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);

h2 = text(13.0,5.1,'TFP  Shocks'); 
h2.Interpreter = 'latex';
h2.FontWeight = 'bold';
h2.FontAngle = 'italic';
h2.FontSize = 20;


subplot(424),plot(dt1,IR2(1:TT1,1+is),trait1,dt1,IR4(1:TT1,1+is),trait2,dt1,IR6(1:TT1,1+is),trait3,'LineWidth',epais),
%axis([0 TT2 -1 4])
title('Oil Inventories','FontWeight','bold','FontSize',mysize)
%xlabel(myxlabel,'FontSize',mysize3,'FontWeight','bold')
ylabel(myylabelshort,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);

subplot(426),plot(dt1,IR2(1:TT1,1+ipo),trait1,dt1,IR4(1:TT1,1+ipo),trait2,dt1,IR6(1:TT1,1+ipo),trait3,'LineWidth',epais),
title('Oil Prices','FontWeight','bold','FontSize',mysize)
%axis([0 TT2 -55 -30])
%xlabel(myxlabel,'FontSize',mysize,'FontWeight','bold')
ylabel(myylabelshort,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);

subplot(428),plot(dt1,IR2(1:TT1,1+ipo)-IR2(1:TT1,1+ipof4),trait1,dt1,IR4(1:TT1,1+ipo)-IR4(1:TT1,1+ipof4),trait2,dt1,IR6(1:TT1,1+ipo)-IR6(1:TT1,1+ipof4),trait3,'LineWidth',epais),
title('Spot minus Futures','FontWeight','bold','FontSize',mysize)
%axis([0 TT2 -0.5 1.0])
xlabel(myxlabel,'FontSize',mysize,'FontWeight','bold')
ylabel(myylabelshort,'FontWeight','bold','FontSize',mysize)
set(gca,'FontSize',mysize);
xt = get(gca, 'YTick');
xt1 = fortick(xt,fmt);
set(gca, 'YTickLabel', xt1, 'YTick', xt);




h = legend(legende1,legende2,legende3);
h.Interpreter = 'latex';
h.FontWeight = 'bold';
h.FontSize = 20;
h.Orientation = 'horizontal';







